% \documentclass[aip,jcp,preprint,unsortedaddress,a4paper,onecolum]{revtex4-1}
\documentclass[aip,jcp,a4paper,reprint,onecolumn]{revtex4-1}
% \documentclass[aps,pre,twocolumn]{revtex4-1}
% \documentclass[aps,jcp,groupedaddress,twocolumn,unsortedaddress]{revtex4}

\usepackage[fleqn]{amsmath}
\usepackage{amssymb}
\usepackage[dvips]{graphicx}
\usepackage{color}
\usepackage{tabularx}
\usepackage{algorithm}
\usepackage{algorithmic}

\makeatletter
\makeatother

\newcommand{\recheck}[1]{{\color{red} #1}}
\newcommand{\redc}[1]{{\color{red} #1}}
\newcommand{\bluec}[1]{{\color{blue} #1}}
\newcommand{\greenc}[1]{{\color{green} #1}}
\newcommand{\vect}[1]{\textbf{\textit{#1}}}
\newcommand{\dd}[1]{\textsf{#1}}


\begin{document}

\title{Report: Accelerated molecular dynamics simulation with optimal control}
% \author{Han Wang}
% \affiliation{Institute for Mathematics, Freie Universit\"at Berlin, Germany}
% \author{Carsten Hartmann}
% \affiliation{Institute for Mathematics, Freie Universit\"at Berlin, Germany}
% \author{Christof Sch\"utte}
% \affiliation{Institute for Mathematics, Freie Universit\"at Berlin, Germany}
% \author{Luigi Delle Site}
% \email{luigi.dellesite@fu-berlin.de}
% \affiliation{Institute for Mathematics, Freie Universit\"at Berlin, Germany}

\begin{abstract}
\end{abstract}

\maketitle

\noindent
To save time, in this report,
I only shortly describe the equations of the optimal control problem, and
present the simulation result for butane.\\

\section{Preliminaries}
\noindent
We follow the work from Hartmann and
Sch\"utte~\cite{hartmann2012efficient}, and assume the Brownian
dynamics of the system:
\begin{align}
  \label{eq:1}
  \dd d\vect x_t =
  -\frac{\sigma^2}{2}\nabla V(\vect x_t)\, \dd dt +
  \beta^{-1/2}\sigma\,\dd d \vect w_t
\end{align}
For simplicity, a homogeneous and scalar $\sigma$ is assumed.
The observable is denoted by $f(\vect x_t)$, and the first mean passage time (FMPT) of some metastable state is denoted
by $\tau$, then the ``work'' of the observable is defined by
\begin{align}
  \label{eq:2}
  W = \int_0^\tau f(\vect x_t)\,\dd dt
\end{align}
The cumulant generating function is given by
\begin{align}
  \label{eq:3}
  J(\vect x) = -\frac1\beta \log \mathbb E_{\vect x}^P [e ^{-\beta W}]
\end{align}
The superscript $P$ means the conditional expectation $\mathbb E_{\vect x}$ is calculated in equilibrium.
It can be shown that the generating function $J$ can be calculated by solving a stochastic optimal
control problem:
\begin{align}
  \label{eq:4}
  J(\vect x) = \min_{\vect u\in U} \mathbb E_{\vect x}^Q
  \big[
  \int_0^\tau f(\vect x_t) + \frac12 \vert \vect u_t\vert^2\, \dd dt\,
  \big]
\end{align}
where the stochastic process $\vect x_t$ is generated by the controlled Brownian dynamics:
\begin{align}
  \label{eq:5}
  \dd d\vect x_t =[ -\frac{\sigma^2}{2}\nabla V(\vect x_t)\, + \sigma \vect u_t\,]\, \dd dt + \beta^{-1/2}\sigma\,\dd d \vect w_t,
\end{align}
and the expectation $ \mathbb E_{\vect x}^Q$ is measured with respect to the trajectory ensemble
generated by Eq.~\eqref{eq:5}.
The running cost $f(\vect x) + \frac12 \vert \vect u\vert^2$ in Eq.~\eqref{eq:4}
is usually denoted by $g(\vect x, \vect u)$.
It can be further proved that the optimal control that minimizes the right-hand-side
of Eq.~\eqref{eq:4} has a feedback form:
\begin{align}
  \label{eq:6}
  \vect u^\ast_t = -\sigma\nabla J(\vect x_t)
\end{align}
which means:
\begin{align}
  \label{eq:7}
  J(\vect x)
  &= \mathbb E_{\vect x}^Q
  \big[
  \int_0^\tau f(\vect x_t) + \frac{\sigma^2}2 \vert \nabla J(\vect x_t)\vert^2\, \dd dt\,
  \big]\\ \label{eq:8}
  \dd d\vect x_t
  &= \frac{\sigma^2}{2}[ - \nabla V(\vect x_t)\, - 2 \nabla J(\vect x_t)] \,\dd dt +
  \beta^{-1/2}\sigma\,\dd d \vect w_t,  
\end{align}
Therefore, out of this self-consistent form, designed a iterative scheme:
\begin{align}
  \label{eq:9}
  J^{(i+1)}(\vect x)
  &= \mathbb E_{\vect x}^Q
  \big[
  \int_0^\tau f(\vect x_t) + \frac{\sigma^2}2 \vert \nabla J^{(i)}(\vect x_t)\vert^2\, \dd dt\,
  \big]\\ \label{eq:10}
  \dd d\vect x_t
  &= \frac{\sigma^2}{2}[ - \nabla V(\vect x_t)\, - 2 \nabla J^{(i)}(\vect x_t)] \,\dd dt +
  \beta^{-1/2}\sigma\,\dd d \vect w_t.  
\end{align}
By letting $f(\vect x_t) = \epsilon$, one can easily show that
\begin{align}
  \label{eq:11}
  \frac{\dd d J_\epsilon}{\dd d \epsilon} \Big\vert_{\epsilon=0}
  =
  \mathbb E_{\vect x}^P[\tau]
\end{align}
By small enough $\epsilon$, one has the finite difference form:
\begin{align}
  \label{eq:12}
  \mathbb E_{\vect x}^P[\tau] \approx J_\epsilon / \epsilon
\end{align}\\

\section{Results}

\subsection{Dumbbell molecule in vacuum}
\begin{figure}
  \centering
  \includegraphics[width=0.49\textwidth]{figs.new/dumbbell/fig-noRed-v.eps}
  \includegraphics[width=0.49\textwidth]{figs.new/dumbbell/fig-freez-v.eps}
  \caption{Dumbbell in vacuum: The original bonded potential and the one calculated by
  the optimal control.}
  \label{fig:0.0}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=0.49\textwidth]{figs.new/dumbbell/fig-noRed-j.eps}
  \includegraphics[width=0.49\textwidth]{figs.new/dumbbell/fig-freez-j.eps}
  \caption{Dumbbell in vacuum:
    The cumulant generating function calculated by a
    brutal force equilibrium simulation and by the optimal control
    simulation. }
  \label{fig:0.1}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=0.49\textwidth]{figs.new/dumbbell/fig-freez-old-j.eps}
  \includegraphics[width=0.49\textwidth]{figs.new/dumbbell/fig-freez-j.eps}
  \caption{Dumbbell in vacuum:
    The cumulant generating function calculated by a
    brutal force equilibrium simulation and by the optimal control
    simulation. Left: $\Delta x = 0.02$ nm, right: $\Delta x = 0.01$ nm.}
  \label{fig:0.1}
\end{figure}


\noindent
To verify the algorithm and numerical implementation, we firstly test a
simplest example: a dumbbell molecule connected by a boned potential
shown in Fig.~\ref{fig:0.0}, which is a symmetric double well
potential with respect to the bond length.
In this system, there does not exist any non-bonded interaction.
The equilibrium population
distribution on the bond-length space, however, is not symmetric,
because of the effect of entropy.  Obviously, the configuration with
longer bond length is more favorable than that with shorter bond
length.  The bond length is chosen as the reaction coordinate, and the
control is applied along the bond length to drive the system towards the
configuration with longer bond length.  \\

\noindent
The bond interaction under
the optimal control is also plotted by the solid blue line in Fig.~\ref{fig:0.0}.
The cumulant generating function calculated by the brutal force equilibrium
simulation  and the optimal control are given in Fig.~\ref{fig:0.1}.
Due to the statistical error,
the output of the iterative scheme Eq.~\eqref{eq:9} and \eqref{eq:10}
is oscillatory. The ``min'' and ``max'' of the oscillatory profiles of the
generating functions are plotted in Fig.~\ref{fig:0.1} by the green and
pink lines. One typical generating is plotted by the blue line with statistical
error indicated by the error bar (with \%95 confidence level).
It is not easy to say the optimal control is perfectly achieved.

\section{Butane in vacuum}

\begin{figure}
  \centering
  \includegraphics[width=0.49\textwidth]{figs.butane.vacuum/fig-j.eps}
  \includegraphics[width=0.49\textwidth]{figs.butane.vacuum.freez/fig-j.eps}
  \caption{Butane in vacuum:
    The cumulant generating function calculated by a
    brutal force equilibrium simulation and by the optimal control
    simulation.
    Left: no constrait. Right: forbids rotation.
  }
  \label{fig:1.1}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=0.49\textwidth]{figs.butane.vacuum/fig-v.eps}
  \includegraphics[width=0.49\textwidth]{figs.butane.vacuum.freez/fig-v.eps}  
  \caption{Butane in vacuum:
    The original dihedral potential and the one calculated by
    the optimal control.
    Left: no constrait. Right: forbids rotation.
  }
  \label{fig:1.0}
\end{figure}

\noindent
The butane molecule in vacuum is studied.
In this system, there does not exist any non-bonded interaction.
Three types bonded interactions are modeled, i.e. bond, angle
and dihedral interactions.
The bond and angle interactions are designed to model the
high-frequency vibration of the covalent bonds and angles between bonds.
The dihedral potential is plotted in
Fig.~\ref{fig:1.0}. In equilibrium, the configurations are equally
distributed over configurations $\phi= -60^\circ$, $\phi= 60^\circ$ and
$\phi= 180^\circ$. The transition among these ``states'' are very slow.
Therefore, in this system,
the dihedral angle is considered to be the slow variable,
and is assumed to be the reaction coordinate. 
We want to design the optimal control in the dihedral angle space 
to drive the system to the
trans-butane configuration ($\phi = 180^\circ$).\\

\noindent
The dihedral interaction plus the optimal control is shown by the
solid blue line in Fig.~\ref{fig:1.0}. 
The cumulant generating function is plotted in Fig.~\ref{fig:1.1}.
The red is the brutal force result and the blue is the optimal control result.


\section{Butane in water}
\begin{figure}
  \centering
  \includegraphics[]{figs/fig-j.eps}
  \caption{Comparison of the cumulant generating function calculated
    by the brutal force equilibrium simulation, and the optimal
    controlled simulation.}
  \label{fig:0}
\end{figure}

\noindent
Now come to the butane in SPC/E water system. The force parameters
are the same as the vacuum butane system.
The reaction coordinate is chosen to be the
dihedral angle (here denoted by $\phi$), and the control is also applied on the
dihedral angle. We calculated the FMPT from any starting
angle to the ``trans''-butane configuration ($\phi = 180^\circ$).
The cumulant generating function is plotted in Fig.~\ref{fig:0}.
The parameter $\epsilon$ is set to be 0.5 and 1.0. We see some systematic error
of the optimal control approach. 
The FMPT calculated by the optimal controlled
dynamics is presented in Fig.~\ref{fig:1} left. The FMPT of the optimal
controlled dynamics itself is presented in Fig.~\ref{fig:1} right.
Roughly, one order of magnitude acceleration is achieved.

\begin{figure}
  \centering
  \includegraphics[width=0.48\textwidth]{figs/fig-fmpt.eps}
  \includegraphics[width=0.48\textwidth]{figs/fig-ctr-fmpt.eps}
  \caption{Left: Comparison of the FMPT calculated by the brutal force equilibrium simulation, and the optimal controlled simulation. Right: The FMPT of  the optimal controlled stochastic process.}
  \label{fig:1}
\end{figure}

% \begin{figure}
%   \centering
%   \includegraphics[]{figs/fig-ctr-fmpt.eps}
%   \caption{The FMPT of  the optimal controlled stochastic process.}
%   \label{fig:2}
% \end{figure}



\section*{References}

\bibliography{ref}{}
\bibliographystyle{unsrt}


\end{document}